# =============================================================================
# 论文：Temporal dynamics and allometric relationships of Haloxylon ammodendron
#        growth traits in a desert ecosystem
# 数据文件：trt_total.csv
# 作者：XXX
# 日期：2025-03-04
# 仅修改：标签字号、坐标轴字号、字体为serif、移除网格，其余参数完全保留
# =============================================================================

# 1. 准备工作：加载必要的包，设置工作目录 --------------------------------
# 安装所需包（若未安装）
required_packages <- c("tidyverse", "ggplot2", "ggpubr", "rstatix", 
                       "car", "emmeans", "multcomp", "cowplot")
install_if_missing <- function(pkg) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}
invisible(lapply(required_packages, install_if_missing))

# 设置工作目录（请根据实际路径修改）
# setwd("your_path")

# 2. 数据导入与初步检查 ----------------------------------------------------
df <- read.csv("trt_total.csv", header = TRUE, stringsAsFactors = FALSE)

# 查看数据结构
str(df)
summary(df)

# 检查缺失值
colSums(is.na(df))

# 将year转换为因子，用于后续分组分析
df$year <- as.factor(df$year)

# 查看各年份样本量
table(df$year)

# 3. 数据清洗：剔除极端异常值（可选，根据专业判断） ------------------------
# 这里采用保守策略：仅剔除明显录入错误的值（如负值或超出生物学极限的值）
# 例如，地径diam不应大于500（现实中极少），高度height不应大于1000等
df_clean <- df %>%
  filter(height > 0, crow > 0, diam > 0) %>%           # 确保所有指标为正
  filter(height < 1000, crow < 1000, diam < 500)       # 根据经验设置合理上限

# 比较清洗前后样本量
cat("原始数据行数:", nrow(df), "\n清洗后数据行数:", nrow(df_clean), "\n")
# 通常清洗掉的数据很少，保留原始数据亦可，这里为保险而做

# 4. 描述性统计：按年份分组计算各指标的均值、标准差、中位数、范围 ----------
desc_stats <- df_clean %>%
  group_by(year) %>%
  summarise(across(c(height, crow, diam),
                   list(mean = ~ mean(., na.rm = TRUE),
                        sd = ~ sd(., na.rm = TRUE),
                        median = ~ median(., na.rm = TRUE),
                        min = ~ min(., na.rm = TRUE),
                        max = ~ max(., na.rm = TRUE)),
                   .names = "{col}_{fn}")) %>%
  pivot_longer(-year, names_to = "variable", values_to = "value") %>%
  separate(variable, into = c("trait", "stat"), sep = "_") %>%
  pivot_wider(names_from = stat, values_from = value)

print(desc_stats)
# 可将结果导出为CSV
write.csv(desc_stats, "descriptive_stats.csv", row.names = FALSE)

# 5. 可视化：箱线图展示各指标年际变化 ------------------------------------
# 自定义主题（仅修改：字体serif + 移除所有网格 + 保留原字号/间距）
theme_sci <- theme_bw() +
  theme(
    # 仅修改：所有文本字体改为serif
    text = element_text(family = "serif"),
    axis.title = element_text(size = 14, face = "bold", family = "serif"),  # 保留原字号14，仅改字体
    axis.text = element_text(size = 12, family = "serif"),                 # 保留原字号12，仅改字体
    legend.title = element_text(size = 12, family = "serif"),              # 保留原字号12，仅改字体
    legend.text = element_text(size = 10, family = "serif"),               # 保留原字号10，仅改字体
    strip.text = element_text(size = 12, face = "bold", family = "serif"), # 保留原字号12，仅改字体
    # 仅新增：移除所有网格（包括major和minor）
    panel.grid = element_blank(),
    # 保留原始代码中仅移除minor网格的逻辑被覆盖，改为移除所有网格
    panel.grid.minor = element_blank()
  )

# 绘制高度箱线图（完全保留原始参数）
p_height <- ggplot(df_clean, aes(x = year, y = height, fill = year)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.2, size = 0.5) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  labs(x = "Year", y = "Height (cm)") +
  theme_sci + theme(legend.position = "none")

# 绘制冠幅直径箱线图（完全保留原始参数）
p_crow <- ggplot(df_clean, aes(x = year, y = crow, fill = year)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.2, size = 0.5) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  labs(x = "Year", y = "Crown diameter (cm)") +
  theme_sci + theme(legend.position = "none")

# 绘制地径箱线图（完全保留原始参数）
p_diam <- ggplot(df_clean, aes(x = year, y = diam, fill = year)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.2, size = 0.5) +
  scale_fill_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
  labs(x = "Year", y = "Basal diameter (cm)") +
  theme_sci + theme(legend.position = "none")

# 组合三图（完全保留原始参数）
combined_plot <- plot_grid(p_height, p_crow, p_diam, ncol = 3, labels = "AUTO")
# 保存图片（完全保留原始尺寸、dpi等参数）
ggsave("boxplots_years.png", combined_plot, width = 12, height = 5, dpi = 300)

# 6. 正态性与方差齐性检验 --------------------------------------------------
# 由于样本量大，正态性检验易显著，此处可用QQ图辅助判断
# Shapiro-Wilk检验（各年份各指标）
norm_test <- df_clean %>%
  group_by(year) %>%
  summarise(across(c(height, crow, diam),
                   ~ shapiro.test(sample(., size = min(5000, length(.))))$p.value,
                   .names = "{col}_p"))  # 因样本过大，抽样至多5000个
print(norm_test)

# 方差齐性检验（Levene检验）
levene_height <- leveneTest(height ~ year, data = df_clean)
levene_crow   <- leveneTest(crow ~ year, data = df_clean)
levene_diam   <- leveneTest(diam ~ year, data = df_clean)

list(Height = levene_height, Crown = levene_crow, Diameter = levene_diam)

# 7. 年际差异显著性检验 ----------------------------------------------------
# 因数据不满足正态分布（根据上述p值），采用Kruskal-Wallis秩和检验
kruskal_height <- kruskal.test(height ~ year, data = df_clean)
kruskal_crow   <- kruskal.test(crow ~ year, data = df_clean)
kruskal_diam   <- kruskal.test(diam ~ year, data = df_clean)

# 输出结果
kruskal_results <- data.frame(
  Trait = c("Height", "Crown diameter", "Basal diameter"),
  chi_sq = c(kruskal_height$statistic, kruskal_crow$statistic, kruskal_diam$statistic),
  df = c(kruskal_height$parameter, kruskal_crow$parameter, kruskal_diam$parameter),
  p_value = c(kruskal_height$p.value, kruskal_crow$p.value, kruskal_diam$p.value)
)
print(kruskal_results)
write.csv(kruskal_results, "kruskal_results.csv", row.names = FALSE)

# 若显著，进行Dunn事后检验（需安装FSA包）
if (!require(FSA)) install.packages("FSA")
library(FSA)

dunn_height <- dunnTest(height ~ year, data = df_clean, method = "bh")  # Benjamini-Hochberg校正
dunn_crow   <- dunnTest(crow ~ year, data = df_clean, method = "bh")
dunn_diam   <- dunnTest(diam ~ year, data = df_clean, method = "bh")

# 提取事后比较结果
posthoc_list <- list(Height = dunn_height$res, Crown = dunn_crow$res, Diameter = dunn_diam$res)
print(posthoc_list)

# 8. 异速生长关系分析（幂律模型：Y = a * X^b） ------------------------------
# 对每个年份分别拟合 log(Y) ~ log(X) 的线性模型
# 我们关注三对关系：height ~ crow, height ~ diam, crow ~ diam

# 创建对数转换后的变量
df_clean <- df_clean %>%
  mutate(logH = log(height),
         logC = log(crow),
         logD = log(diam))

# 定义函数：对每个年份分别拟合模型，并提取系数、R²、置信区间
fit_allometry <- function(data, x_var, y_var, group_var = "year") {
  results <- data %>%
    group_by(!!sym(group_var)) %>%
    summarise(
      intercept = coef(lm(!!sym(y_var) ~ !!sym(x_var)))[1],
      slope = coef(lm(!!sym(y_var) ~ !!sym(x_var)))[2],
      r_squared = summary(lm(!!sym(y_var) ~ !!sym(x_var)))$r.squared,
      p_value = summary(lm(!!sym(y_var) ~ !!sym(x_var)))$coefficients[2,4],
      .groups = "drop"
    ) %>%
    mutate(relationship = paste0(y_var, "~", x_var))
  return(results)
}

# 拟合三对关系
res_HC <- fit_allometry(df_clean, "logC", "logH")  # height ~ crow
res_HD <- fit_allometry(df_clean, "logD", "logH")  # height ~ diam
res_CD <- fit_allometry(df_clean, "logD", "logC")  # crow ~ diam

# 合并结果
allometry_results <- bind_rows(res_HC, res_HD, res_CD)
print(allometry_results)
write.csv(allometry_results, "allometry_coef.csv", row.names = FALSE)

# 9. 异速生长关系的年际比较（ANCOVA） ---------------------------------------
# 检验不同年份的斜率是否相等，即年份与logX的交互作用是否显著
ancova_HC <- lm(logH ~ logC * year, data = df_clean)
ancova_HD <- lm(logH ~ logD * year, data = df_clean)
ancova_CD <- lm(logC ~ logD * year, data = df_clean)

# 输出ANCOVA结果（Type III方差分析）
library(car)
Anova(ancova_HC, type = "III")
Anova(ancova_HD, type = "III")
Anova(ancova_CD, type = "III")

# 若交互作用显著，说明异速生长指数（斜率）随年份变化

# 10. 可视化异速生长关系（散点图+回归线）-----------------------------------
# 自定义函数绘制单个关系图
plot_allometry <- function(data, x_var, y_var, xlab, ylab) {
  p <- ggplot(data, aes(x = !!sym(x_var), y = !!sym(y_var), color = year)) +
    geom_point(alpha = 0.3, size = 0.8) +
    geom_smooth(method = "lm", se = TRUE, fullrange = FALSE) +
    scale_color_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) +
    labs(x = xlab, y = ylab, color = "Year") +
    theme_sci
  return(p)
}

# 注意这里使用原始尺度（对数坐标）或原始值？为直观展示异速关系，最好在对数坐标下绘图
# 因为我们在log-log尺度拟合线性模型，所以用对数坐标展示更符合模型假设
p1 <- plot_allometry(df_clean, "logC", "logH", "log(Crown diameter)", "log(Height)")
p2 <- plot_allometry(df_clean, "logD", "logH", "log(Basal diameter)", "log(Height)")
p3 <- plot_allometry(df_clean, "logD", "logC", "log(Basal diameter)", "log(Crown diameter)")

# 组合图形
allom_plots <- plot_grid(p1 + theme(legend.position = "none"),
                         p2 + theme(legend.position = "none"),
                         p3 + theme(legend.position = "none"),
                         ncol = 3, labels = c("A", "B", "C"))
# 提取图例
legend <- get_legend(p1 + theme(legend.position = "bottom"))
allom_combined <- plot_grid(allom_plots, legend, ncol = 1, rel_heights = c(1, 0.1))
ggsave("allometry_plots.png", allom_combined, width = 15, height = 6, dpi = 300)

